Electrostatics in Periodic Boundary Conditions and Real-space Corrections 
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We address periodic-image errors arising from the use of periodic boundary conditions to describe 
systems that do not exhibit full three-dimensional periodicity. The difference between the periodic 
potential, as straightforwardly obtained from a Fourier transform, and the potential satisfying any 
other boundary conditions can be characterized analytically. In light of this observation, we present 
an efficient real-space method to correct periodic-image errors, based on a multigrid solver for the 
potential difi'erence, and demonstrate that exponential convergence of the energy with respect to cell 
size can be achieved in practical calculations. Additionally, we derive rapidly convergent expansions 
for determining the Madelung constants of point-charge assemblies in one, two, and three dimensions. 
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I. INTRODUCTION 

First-principles calculations frequently employ peri- 
odic boundary conditions to predict materials properties. 
Besides constituting a natural choice when studying crys- 
talline systems, periodic boundary conditions allow the 
use of highly optimized fast Fourier transform (FFT) al- 
gorithms [if, [2I, [sj, which considerably reduce the com- 
putational cost associated with the resolution of elec- 
trostatic equations, and allow an efficient evaluation of 
electronic kinetic energies and interatomic forces when 
used in conjunction with a plane- wave basis set. Despite 
these algorithmic advantages, periodic boundary condi- 
tions require large supercells when studying aperiodic or 
partially periodic systems (e.^., isolated molecules, poly- 
mer chains, and slabs) in an effort to minimize spuri- 
ous electrostatic interactions between periodic images Q . 
Charged systems are particularly problematic, since con- 
ventional algorithms automatically enforce charge neu- 
trality by introducing an artificial jellium background 
(Note that the electrostatic energy of a charged system 
exhibiting three-dimensional periodicity is infinite.) As 
shown by Makov and Payne, these artifacts induce sig- 
nificant errors scaling as for the energy of neutral 
polarized systems and 1/L for that of charged systems, 
where L denotes the size of the unit cell [sj]. 

In addition to the Makov-Payne asymptotic correc- 
tion [5[, several schemes have been devised to reduce 
periodic-image errors. Barnett and Landman proposed 
to eliminate periodic-image interactions for cluster sys- 
tems by restricting the plane-wave expansions of the 
wavefunctions and of the charge density to a spherical 
domain in reciprocal space [1, 0, [1] . A generalization of 
this reciprocal-space a^)roach was introduced by Mar- 
tyna and Tuckerman [9]. The electrostatic-cutoff ap- 
proach proposed by Jarvis, White, Godby, and Payne 
suppresses periodic-image effects by damping the elec- 



* Electronic address: dabo'is^mit.edul 



trostatic potential beyond a certain interaction range 
[13]. The corrective method introduced by Blochl con- 
sists of using atom-centered Gaussian charges and Ewald 
summation techniques to cancel periodic-image interac- 
tions [11]. In the local- moment-countercharge (LMCC) 
method developed by Schultz, a superposition of Gaus- 
sians is employed as a local-moment model for calculating 
the Coulomb potential analytically up to a certain multi- 
pole order, the remaining electrostatic contribution being 
computed using conventional plane- wave techniques [12]. 
Considering atomic adsorption on neutral slabs, Neuge- 
bauer and Scheffier proposed eliminating the adsorbate- 
induced polarization through the introduction of a coun- 
teracting planar dipole between slab images [14]. Refine- 
ments of this method, based on the linear- and planar- 
average approximations proposed by Baldereschi, Baroni, 
and Resta [l^ , were subsequently developed [13, E, [l^ • 
Extending this approach to charged surfaces, the pre- 
scription of Lozovoi and Alavi relies on inserting a Gaus- 
sian layer in vacuum to compensate for the excess charge 
and to allow electric-field discontinuities across the layer 

a. 



In this work, we propose an alternative approach for 
correcting periodic-image errors and show that exponen- 
tial energy convergence with respect to cell size can be 
obtained at tractable computational cost. The approach 
proceeds by calculating the electrostatic potential in real 
space, exploiting the periodic solution of the Poisson 
equation computed using inexpensive FFT techniques. 
In the following sections, we first discuss and charac- 
terize the difference between the open-boundary electro- 
static potential and its periodic counterpart, providing a 
comparative basis for analyzing the relative accuracy of 
various corrective schemes. Second, we present our cor- 
rection method and assess its performance. Last, we ex- 
tend the method to the study of systems exhibiting one- 
or two-dimensional periodicity, beyond the conventional 
linear- and planar- average approximations. 
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II. COMPARISON OF THE OPEN-BOUNDARY 
AND PERIODIC POTENTIALS 

A. Definition of the Corrective Potential 

The electrostatic potential v generated by a charge dis- 
tribution p satisfies the Poisson equation: 

V^v(r) = -47rp(r) (1) 

(atomic units are used throughout). In the absence of 
an external electric field, we can solve Eq. [T] subject to 
open-boundary conditions ('i;(r) ^ as |r| ^ +oo). As a 
result, the electrostatic potential v can be computed via 
Coulomb integration: 

(Although this study focuses on open boundary condi- 
tions, it should be noted that the contribution from an 
external field E can be incorporated by adopting the 
asymptotic boundary conditions v{y) — E • r, which 
simply adds a term — E • r to the solution of the Pois- 
son equation.) A differential equation similar to Eq. [1] 
can be written for the periodic potential v' ^ keeping in 
mind that periodic boundary conditions can only accom- 
modate a net zero charge (as seen from Gauss' law): 

V\'{v) = -47r(p(r) - {p)). (3) 

As a consequence, the periodic potential can be evaluated 
in the reciprocal-space representation as: 

= E (4) 

where we set the arbitrary component v' = 0) = {v') 
to zero. 

It should be noted that the open-boundary potential 
V and its periodic counterpart v' are distinct. We define 
the corrective potential v^'^'^^ as the difference v — v'. The 
potential v^'^'^^ must satisfy: 

V2t;--(r) = -47r(p), (5) 

for which we specify Dirichlet boundary conditions: 
^corr _ ^ _ ^YiQ cell boundaries. (Note that the 

solution of this elliptic boundary value problem [l5|, [l^ 
is uniquely defined.) Eq. [5] indicates that the curvature 
of the corrective potential is a constant. It should also be 
noted that, apart from the value of the average (p), Eq. [5] 
is independent of the structural details of the charge den- 
sity p. Instead, these details are entirely embedded in the 
Dirichlet boundary conditions, which reflect the electro- 
static contributions from compensating jellium and from 
the surrounding images. 

In order to illustrate the implications of Eq. [5l we 
consider a pyridazine cation in a periodically repeated 



cubic cell. The open-boundary potential the periodic 
potential v' ^ and the corrective potential ^i^^^^ are shown 
in Figure [TJ First, we observe that the potential v' is 
shifted down in energy with respect to due to the fact 
that the average {v') is null by construction. In addition 
to this energy shift, the potential v' is significantly dis- 
torted. This distortion results from satisfying the period- 
icity conditions. Most importantly, we observe that the 
corrective potential v^^^^ varies smoothly over space. The 
smooth spatial dependence of '^^^^^ contrasts markedly 
with the strong variations in v and in v' . Performing a 
polynomial regression, we can verify that the potential 
^corr -g quadratic to good approximation in the proxim- 
ity of the cell center with departures from parabolicity 
restricted to the vicinity of the periodic boundaries. 

To further examine the characteristics of v^^^^ ^ we con- 
sider the adsorption of carbon monoxide molecules on 
neutral and charged platinum slabs. Following Neuge- 
bauer and Scheffler, the electrostatic correction is cal- 
culated along the z-direction within the planar-average 
approximation (that is, from the x?/- aver age of the charge 
distribution) [l^]. The validity of this approximation is 
discussed in the last section. For CO molecules adsorbed 
on a neutral slab (Figure [2^), the periodic potential is 
shifted up in energy and tilted with respect to the open- 
boundary potential. The potential correction is seen to 
be linear, in agreement with the analysis of Neugebauer 
and Scheffler [ij]. For CO molecules adsorbed on a slab 
of surface charge d (Figure [JJd), the real-space potential 
diverges as 4:7ra\z\. In this case, the periodic potential 
undergoes a significant energy downshift, which decreases 
the energy of the positively charged slab. Moreover, we 
observe that v' is significantly curved in the slab region. 
Consistent with these observations and with Eq. O the 
corrective potential v^^^^ is found to be parabolic every- 
where in the unit cell. 



B. Quasiparabolic Behavior of the Corrective 
Potential 

In order to complete the analysis of the corrective po- 
tential, we consider a point charge g = +e in a peri- 
odically repeated cubic cell of length L, as illustrated 
in Figure [3l The corrective potential generated by the 
uniform jellium and the surrounding point charges is de- 
noted Vq^^^. Note that Vq^^^ cannot be calculated di- 
rectly as the difference between the potential of a lattice 
of point charges Vq and the point-charge potential 1/r 
since the representation of a point charge in reciprocal 
space requires an infinite number of plane-wave compo- 
nents. Instead, to obtain vg^^^, we can exploit the cubic 
symmetry of the system, writing the corrective potential 
as: 

i;g^"^(r) = ^g^""(r = 0) 

+VV^^(r = 0)^ + O(|rr). (6) 
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(a) v(v)_ (b) t/(r) 




FIG. 1: (a) Open-boundary electrostatic potential v, (b) periodic electrostatic potential v' , and (c) electrostatic-potential 
correction y^^^^ — y — y' for a pyridazine cation in a cubic cell of length L = 15 bohr. The potentials are plotted in three 
orthogonal planes (Oxy), (Oxz), and iOyz) passing through the center of the cell. 



This parabolic expansion, valid up to third order, con- 
firms that the point-charge correction Vq^^^ is almost 
quadratic in the vicinity of r = 0. For noncubic lattices, 
due to inversion symmetry, the point-charge corrective 
potential takes a more general form: 




(r = 0)r^ + O(|r|^), (7) 



where (re) are the coordinates of r along the principal 
axes. Thus, the corrective potential in a noncubic lattice 
is also quasiparabolic. 

Turning now to an arbitrary distribution p, we can ex- 
press the electrostatic correction ^^^^^ by superposition: 



^corr^^^ = / ^g^^^(r - Y')p{Y')dY'. (8) 



As a consequence, defining Vmax as the distance beyond 
which the parabolic expansion (Eq. [6]) ceases to be valid, 
the corrective potential ^i^^^^ can be considered as nearly 
parabolic, provided that the spread of the distribution is 
tolerably lower than Vmax- 



C. Connection with Existing Schemes 

Having justified the general characteristics of the 
electrostatic-potential correction, we now determine the 
terms in the expansion of Vq^'^^ (Eq. [6]). The potential 
at the origin VQ^^^(r = 0) can be written in terms of 
the Madelung constant ao [22|] of a cubic lattice of point 
charges in a compensating jellium background: 

^g^"^(r = 0) = ^. (9) 

(The calculation of the Madelung constant of a jellium- 
neutralized assembly of point charges is discussed in Ap- 
pendix [A)) Note that VQ^^^{r = 0) is positive, reflecting 
the stabilizing contribution from the jellium compensa- 
tion. The value of \/'^VQ^'^^{r = 0) is then determined 
from Eq. [5} 

4-7r 

V'v^'-ir = 0) = -j^. (10) 
Hence, the point-charge correction can be expanded as: 

^corr^r) = ^-^r^ + 0{\r\% (H) 

The terms in this parabolic expansion bear a strong re- 
semblance to those entering into the Makov-Payne cor- 
rection [5]. This correspondence is discussed further in 
Sec. [HDl 



4 



-1 - 





10 



7.5 - 



I 2.5 
o 




- 



10 20 30 40 50 

coordinate z (bohr) 

FIG. 2: Open-boundary electrostatic potential v, periodic 
potential v\ and electrostatic-potential correction y^^^^ av- 
eraged in the x^-plane parallel to the surface for (a) carbon 
monoxide adsorbed on a neutral platinum slab, and (b) car- 
bon monoxide adsorbed on a charged platinum slab. 



The above expansion allows us to approximate the elec- 
trostatic correction induced by a set of compensating 
charges. Indeed, introducing N charges, we can define 
a parabolic point-countercharge (PCC) potential Vp^'C 
as: 



PCcir) = T.1n(^-§^,{r-r^) 

n=l ^ 



This expression may be rewritten: 



ooq 
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(12) 



(13) 



where q = ^^Qn is the total charge, p = ^^qn^n 
denotes the total dipole moment, and Q = Qn^n 
stands for the total quadrupole moment of the coun- 
tercharge distribution. Eq. [13] indicates that parabolic 
PCC schemes can correct periodic-image errors up to 
quadrupole- moment order. Note that no more than 
Nmax = 7 countercharges are sufficient to obtain the 
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FIG. 3: Corrective potential Vq^^^ for a cubic lattice of point 
charges and its parabolic approximation in the vicinity of the 
origin. 



most accurate parabolic correction (one charge for q^ two 
for p, and four for Q). To obtain higher-order PCC 
corrections, one would need to determine more terms in 
the expansion of the point-charge correction, beyond the 
parabolic contributions. An example of accurate calcu- 
lations using harmonic expansions can be found in Ref. 

m- 

An alternative approach is to employ countercharges 
whose corrective potential can be computed handily. A 
popular choice is to use Gaussian densities, as proposed 
by Blochl [ll|. Repeating the preceding analysis for a 
Gaussian density of charge q = -\-e^ we can expand the 
Gaussian corrective potential v^^l^ as: 



27r 
31? 



r'^0{\r\% 



(14) 



where acr/L is the Madelung constant of an assembly of 
Gaussians of width <j immersed in a compensating jellium 
in a cubic cell of length L. It is more convenient, however, 
to write the corrective potential directly as: 



(r) 



= '-^^ - V ije-^^VVs.r (15) 

where Vcr is the electrostatic potential of an isolated Gaus- 
sian charge, and v'^ ^ is the potential corresponding to a 
periodically repeated Gaussian in a jellium background. 
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FIG. 4: Point- countercharge (PCC), Gaussian-countercharge (GCC), and density- countercharge (DCC) corrective potentials 
for a pyridazine cation C4H5NJ in a cubic cell of length L = 15 bohr. The corrective potentials are plotted along the z-axis 
perpendicular to the plane of the molecule, as defined in Figured] The PCC and GCC corrections are calculated up to dipole 
order. The spread of the Gaussian countercharges is a = 0.5 bohr. 



The sum in the right-hand side of the equation con- 
verges very rapidly, and can be calculated using FFT 
techniques. Superimposing N compensating charges, 
the Gaussian-countercharge (GCC) corrective potential 
^GCC expressed as: 

N 

«GCc(r) = E9»<rr(r-r„). (16) 

n=l 

This results in the following approximation for the open- 
boundary potential v: 

v{r)^v'{r)+v^accir)- (17) 

We underscore that this scheme is equivalent to the Gaus- 
sian scheme introduced by Blochl [11] and the LMCC 
method proposed by Schultz [12]. The equivalence with 
LMCC approach can be established by recasting Eq. [T71 
as: 

v{t) ^ vpBc(r) + VGCc(r) 

(18) 

vpBc{r) =v'{r) -v'^cci^)^ 

where vgcc{'^) = ^Qn^ai'^ — ^n) is the electrostatic po- 
tential generated by the isolated countercharge distribu- 
tion, and VQQQ{r) = ^ ^n'^^ ~ ^n) is the correspond- 
ing periodic potential. 

We are now in a position to compare the correc- 
tive potentials '^p^'c ^gcc ^i^^ potential v^^^^^ 



obtained as the direct difference between the open- 
boundary potential and its periodic counterpart. For 
our comparative analysis, we refer to the exact correc- 
tive potential v^^^^ as the density-countercharge (DCC) 
potential. The DCC potential is obtained by evaluating 
the Coulomb integral defining v at each grid point in the 
unit cell. (A cheaper alternative to this procedure is pre- 
sented in the next section.) The PCC, GCC, and DCC 
potentials for a charged pyridazine cation in a cubic cell 
of length L = 15 bohr are plotted in Figure [H The PCC 
and GCC corrections are computed up to dipole order. 
First, it should be noted that the maximal energy of the 
PCC potential is slightly above its GCC counterpart, re- 
flecting the fact that the Madelung energy of an array of 
point charges immersed in a jellium is higher than that 
of a jellium-neutralized array of Gaussian charges (cf. 
Appendix [A|) . In addition, the maximal DCC energy is 
found to be approximately 0.05 Ry above aoq/L, indi- 
cating that the dipole PCC and GCC corrections tend 
to underestimate the energy of the system. Moreover, 
the parabolic PCC potential is not as steep as its GCC 
counterpart, suggesting that the energy underestimation 
will be more significant for the GCC correction. Owing 
to the cubic symmetry of the cell, the PCC and GCC 
potentials display the same curvature in each direction 
of space, equal to one third of — 47r(p). In contrast, the 
curvature of the DCC potential is not uniform, due to 
the nonspherical nature of the molecular charge density. 
This shape dependence suggests that the accuracy of the 
GCC correction could be improved by optimizing the ge- 
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FIG. 5: Electrostatic energy of two Gaussians of unit charge 
and unit spread calculated via (a) real-space integration, (b) 
reciprocal-space integration with the PCC energy correction 
given by Eq. [21] and (c) reciprocal-space integration with the 
energy correction given by Eq. 15 in Ref. [5]. The Gaussian 
charges are positioned at ro = (—5, —5, —5) and ri = (5, 5, 5) 
(corresponding to a quadrupole moment Q of 153 a.u.). 

ometry of the Gaussian countercharges. 

In summary, we have shown that the PCC (Makov- 
Payne), GCC (LMCC), and DCC corrections belong 
to the same class of periodic-image corrections. The 
analysis of the corrective potential has established that 
the parabolic PCC correction cannot eliminate periodic- 
image interactions beyond quadrupole order. Difficulties 
inherent in the GCC scheme have also been evidenced. 
To overcome these limitations, an efficient implementa- 
tion of the DCC correction is presented in Sec. Illli 
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FIG. 6: DCC total energy as a function of cell size for a 
pyridazine cation varying the coarse-grid cutoff E^uV from 
10 Ry (M = 13 x 13 x 13) to 250 Ry (M = AT = 73 x 73 x 73). 
Also depicted is the corrective potential y^^^^ in the plane of 
the molecule as a function of the coarse-grid resolution at a 
cell size of 15 bohr. 



D. Energy Correction 

To conclude this preliminary analysis, we give the ex- 
pression of the energy correction A£^^^^^ in terms of the 
corrective potential v'^'^^^. The total electrostatic energy 
of the system being equal to: 

E=\J vir)pir)dr, (19) 

the corrective energy can be expressed as fv^ : 

A^^^^^ = ^J ^^^^^(r)p(r)dr. (20) 

it is worth mentioning that in the case of a single point 
countercharge Q = J p(r)(ir, the PCC energy correction 
can be written as: 

^^eorr ^ ^ qv^'^^^ [v) p{v)dv 

~ 2L 3L3- ^^^> 



The first term corresponds to the Madelung energy cor- 
rection, as proposed by Leslie and Gillian [4]. Note that 
the second term differs from Eq. 15 in Ref. ^] by a factor 
1/2. The validity of the energy correction given by Eq. 
[21] is illustrated in Figure [5l 



III. IMPLEMENTATION OF THE 
DENSITY-COUNTERCHARGE CORRECTION 

A. Density-countercharge Algorithm 

In the preceding section, the corrective potential 
ycorr _ ycorr^ ^as Calculated directly by subtracting 
the periodic potential from its open-boundary counter- 
part. The computational cost of this direct method is 
prohibitively high, on the order 0{N'^) (where TV is the 
number of grid points), corresponding to the evaluation 
of Coulomb integrals at each point of the grid. In this 
section, we present a scheme that reduces this computa- 
tional burden. The scheme exploits both the Poisson 
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equation for v^^^^ (Eq. [5]) and the fact that v^^^^ is 
smoothly varying. 

First, we note that taking into account appropriate 
boundary conditions, Eq. [5] can be solved efficiently us- 
ing multigrid solvers [H, Q, [H, Hi, [13, H^. Multigrid 
algorithms typically scale as 0{N log N), that is, com- 
parable to the scaling of an FFT computation. Hence, 
the overall cost of the calculation can be reduced from 
0{N'^) to 0(A^^/^), corresponding to the expense arising 
from the determination of the boundary conditions. Al- 
though a similar approach may be employed to directly 
solve the electrostatic equation defining v (Eq. [1]), we 
emphasize that Eq. [5] allows a considerable reduction in 
numerical error in the finite-difference evaluation of the 
electronic Laplacian — since ^^^^^ is much smoother than 

V. 

Further exploiting this idea, it is possible to solve Eq. 
[5] on a grid much coarser than that used to discretize 
the charge density. To illustrate this fact, we consider a 
pyridazine cation in a periodic cubic cell of varied size 
(Figure [6]). The total energy of the system is calculated 
using density- functional theory [l^. An energy cutoff 
Ecut = 250 Ry is applied to the plane-wave expansion of 
the charge density. The total energies are corrected using 
the DCC scheme by solving the electrostatic equation of 
ycorr ^ coarse grid for several values of the energy 
cutoff, denoted E^^^^. Reducing the energy cutoff E^'^l^ 
from 250 to 40 Ry, the corrected energies are observed 
to depart by less than 5 x 10~^ Ry from their converged 
values for cell sizes greater than 13 bohr. The ability to 
decrease the number of grid points without a significant 
loss of accuracy enables a substantial reduction of the 
additional computational cost from 0{N^^^) to 0(M^/^), 
where M is the number of coarse-grid points. Note that 
diminishing the plane-wave energy cutoff from 250 to 40 
Ry at L = 15 bohr reduces the cost of the boundary- 
condition calculation by a factor 29^/73^ ^ 1/100. 

Before presenting the algorithm, we draw attention to 
the fact that the DCC scheme relies on the central idea 
that most of the structural characteristics of the open- 
boundary potential v can be removed by subtracting out 
its periodic counterpart v' . The residual v'^^^^ (that is, 
the amount by which v' fails to reproduce v) is smooth 
and can be determined on a coarse grid at low computa- 
tional cost. Additional computational savings come from 
the ability to avoid updating the potential v^^^^ at each 
step of the self-consistent-field (SCF) calculation, but in- 
stead at fixed interval between electronic iterations. 

The DCC algorithm for a typical electronic-structure 
calculation can be described as follows. Let 7V^^^^ denote 
the number of SCF steps between each update of the 
corrective potential. 

1. Start from an initial charge distribution p on the 
fine grid. 

2. Calculate the periodic potential corresponding 
to p. 



3. Transfer p and on the coarse grid (tricubic in- 
terpolation [3Q,]) to obtain the coarse-grid density 
p and coarse-grid periodic potential v' . 

4. Calculate the real-space potential v at the bound- 
aries of the coarse grid from p to obtain the Dirich- 
let boundary conditions v^^^^ = v — v' . 

5. Solve y'^qj^orr _ _4^^p^ (multigrid techniques) to 
obtain the corrective potential i)^^^^. 

6. Transfer i;^^^^ on the fine grid (tricubic interpola- 
tion) to obtain '^^^^^^ and calculate v = v^'^'^'^ + v' . 

7. Perform 7V^^^^ electronic SCF steps. 

8. Iterate from Step 2 until reaching SCF convergence. 

Note that we employ real-space tricubic interpolation 
techniques in order to avoid oscillatory distortions in- 
herent in Fourier- transform interpolation schemes. We 
also underscore that the DCC algorithm can be efficiently 
parallelized, since its most expensive step (namely, the 
calculation of the Dirichlet boundary conditions) scales 
linearly with the number of processors. 

The above procedure can be adapted to one- and two- 
dimensional systems by considering the linear or planar 
average of the charge density for calculating the correc- 
tive potential [l^. (The validity the linear- or planar- 
average approximations will be discussed in the final sec- 
tion.) The computational cost of this approach is mod- 
erate, on the order of 0(M^/^) and 0{M) for one and 
two dimensions, respectively. 

It should also be mentioned that the DCC algorithm 
can be used in combination with multipole-expansion 
methods for a rapid evaluation of the Dirichlet bound- 
ary conditions (Step 4). The accuracy of this approach 
depends on the precision of the multipole expansion at 
the boundary of the supercell. (A mathematical discus- 
sion on the long-range accuracy of multipole expansions 
is presented in Sec. 3.4. of Greengard's dissertation [^.) 
The performance the multipole-expansion approach is re- 
ported in Appendix [Bl 



B. Applications 

The energy of a pyridazine molecule as a function of 
cell size L for each countercharge correction is reported in 
Figure [71 For this neutral species, the uncorrected energy 
shows a characteristic minimum at L = 14 bohr before 
slowly approaching its asymptotic value. In contrast, the 
corrected energies are seen to converge monotonically to- 
wards their common energy limit. Although the three 
schemes demonstrate comparable convergence, it should 
be noted that the PCC method is slightly more accurate. 
In addition to further validating the energy expansion 
given by Eq. [211 this comparison suggests that the PCC 
correction can be preferred for studying neutral species. 
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FIG. 7: Total energy as a function of cell size for a neutral 
pyridazine molecule without correction and corrected using 
the PCC, GCC, and DCC schemes. The PCC and GCC cor- 
rections are calculated up to quadrupole order. The inset 
shows a pyridazine molecule in a cell of size L = 15 bohr. 



with the notable exception of elongated systems (e.^., 
polymer fragments or terminated nanotubes). 

We now consider the energy of a pyridazine cation as a 
function of cell size (Figure [8]). We use energy cutoffs of 
35 and 250 Ry for expanding the wavefunctions and the 
charge density, and select a coarse-grid cutoff of 35 Ry for 
calculating the DCC correction. Expectedly, the uncor- 
rected energy converges very slowly with respect to L (at 
19 bohr, the energy error is still larger than 0.15 Ry). The 
PCC and GCC corrections substantially improve the con- 
vergence of the total energy, reducing periodic-image er- 
rors by one order of magnitude. Using the DCC scheme, 
the energy is observed to converge even more rapidly, re- 
flecting the exponential disappearance of energy errors 
arising from charge density spilling across periodic cells: 
at a cell size of 15 bohr, which is barely larger than the 
size of the molecule, the DCC energy is converged within 
10~^ Ry. The performance of each scheme as a function 
of the total computational time is shown on a logarithmic 
energy scale in Figure [9l Each curve corresponds to cell 
sizes in the range 12-19 bohr. For meaningful compari- 
son with the DCC scheme, the PCC and GCC corrective 
potentials are also updated at fixed SCF intervals. We 
observe that the computational cost of the corrected cal- 
culations is comparable to that without correction for a 
considerable improvement in accuracy. For this charged 
system, the DCC approach constitutes the most advanta- 
geous alternative, improving the energy precision by two 
orders of magnitude over the PCC and GCC corrections 
for cell sizes above 15 bohr. 

The performance of the DCC and GCC corrective 
schemes for a neutral polyvinylidene fluoride (PVDF) 
chain is reported in Figure [TOl The comparison shows a 
significant improvement of energy convergence for both 
schemes. As shown in the inset, the performance of 
the DCC scheme is perceptibly superior to that of the 
GCC scheme. We emphasize that for systems exhibit- 
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FIG. 8: Total energy as a function of cell size for a pyridazine 
cation without correction and corrected using the PCC, GCC, 
and DCC schemes. The PCC and GCG corrections are cal- 
culated up to quadrupole order. The inset of the top graph 
shows a pyridazine cation in a cell of size L = 15 bohr. 



ing one dimensional periodicity, the additional computa- 
tional cost due to the electrostatic correction is moderate, 
on the order of 0(M) at most. 

The DCC scheme can also be used in the calculation 
of work functions, as it solves energy-reference issues by 
automatically setting the vacuum level to zero. Figure 
[TTl depicts the convergence of the opposite Fermi energy 
of a Ft (100) slab as a function of transverse cell size. 
The wavefunction, charge-density, and corrective poten- 
tial energy cutoffs are 25, 200, and 150 Ry, respectively. 
We use a shifted 5x5x1 mesh with a cold-smearing oc- 
cupation function [sg'] (smearing temperature of 0.03 Ry) 
to sample the Brillouin zone. Without correction, the rel- 
ative error in the Fermi energy stays above 100% for all 
cell sizes in the considered range. Using the GCC scheme. 
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1000 1500 
CPU time (s) 

FIG. 9: Accuracy of the total energy as a function of compu- 
tational time without correction and using the PCC, GCC, 
and DCC schemes for cell sizes in the range 12-19 bohr. For 
each scheme, the corrective potential is updated every five 
SCF iterations. 



the convergence of the Fermi level improves greatly: at 
150 bohr, the relative error reduces to approximately 0.1 
eV. Employing the DCC corrective scheme, the calcu- 
lated Fermi energy is converged within 2 meV at 60 bohr 
and 0.1 meV at 150 bohr. Thus, the DCC scheme al- 
lows to directly determine the work function of a metal 
as the opposite of the calculated Fermi energy using su- 
per cells of minimal size. A similar convergence improve- 
ment is obtained for the work function of carbon nan- 
otubes [35[. Besides improving the convergence of total 
energies, the DCC approach can be employed to correct 
structural and vibrational properties [32], and to calcu- 
late linear-response characteristics with a reduced com- 
putational effort [3i,[33,[33. 



IV. BEYOND THE LINEAR- AND 
PLANAR- AVERAGE APPROXIMATIONS 

A. Treating Systems with Partial Periodicity 



-122.64 




10 12.5 15 17.5 20 22.5 
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FIG. 10: Total energy as a function of transverse cell size for 
a polyvinylidene fluoride (PVDF) chain without correction, 
and using the GCC and DCC schemes. 
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FIG. 11: Convergence of the opposite Fermi energy —cf as 
a function of transverse cell size for a Pt(lOO) slab without 
correction, and using the GCC and DCC schemes. 



In the preceding sections, we have assumed that the 
corrective potential of a one- or two-dimensional system 
can be obtained by homogenizing the system along its pe- 
riodicity directions, as initially proposed by Baldereschi, 
Baroni, and Resta [13]. This approach, referred to as 
the linear- or planar- average approximation, has been 
frequently emp loyed in electronic-structure calculations 

0, El, [13, Ska]. 

Alternative schemes adapting the Ewald method to 
evaluate conditionally convergent lattice sums [38] or 
generalizing the FMM approach [39*, ^40] have also been 
proposed for systems exhibiting partial periodicity. Such 
schemes are particularly suited to localized-orbital cal- 
culations but are of relatively limited applicability for 
plane-wave implementations. Here, we propose an effi- 
cient method to calculate the electrostatic potential for 
partially periodic systems, taking into account the full 
three-dimensional structure of the charge distribution. 
In addition to presenting this methodological extension, 
we discuss how to assess the validity of the linear- and 
planar- average approximations a priori in terms of struc- 
tural characteristics of the system. 

B. DCC Scheme for One-dimensional Periodicity 

To introduce the DCC approach for one-dimensional 
systems, we first study the electrostatic problem corre- 
sponding to an isolated sinusoidal-density line: 

pir) = S^^\r^)eMWzz), (22) 

where J^^^ stands for the two-dimensional Dirac delta 
function and denotes the transverse coordinates (x, y). 
Making the ansatz v{r) = Q{r±;gz)ex.-p{igzz) for the 
Green's function, we obtain: 

{Vl-gl)g{rr,9.) = -47r5(2)(r^). (23) 
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FIG. 12: Fourier-decomposition calculation of the electrostatic potential ^;(r^, z) = v{r±; gz)e^^^^ for an infinite polyvinyli- 
dene fluoride (PVDF) chain. (1) The longitudinal Fourier transform of the charge density is calculated to obtain the contribu- 
tions from each axial wavevector Qz; (2) the electrostatic potential generated by each Fourier component of the charge density 
is calculated using Green's functions; (3) the electrostatic potential is then transformed back to real space. 



The solution of this generalized electrostatic problem can 
be written as: 

r a(rx;0) = -21n|rx|, 

{ (24) 
[ G{rr.g,) = 2Ko{g,\r^\) for ^ 

where Kq is the modified Bessel function of the second 
kind. Note that Ko{gz\T±\) = — ln|rx| + ... when ^^Ir^l 
approaches zero, reflecting the fact that a sinusoidal- 
density line can be considered as uniform when seen from 
a distance much smaller than its wavelength. Know- 
ing the electrostatic potential generated by a single line 



(the Green's function characterizing the generalized elec- 
trostatic problem), the potential of an arbitrary one- 
dimensional charge distribution can be determined an- 
alytically, as illustrated in Figure [121 The general proce- 
dure consists of calculating the one-dimensional Fourier 
transform of p to obtain its longitudinal Fourier compo- 
nents p{r_\_;gz) (step 1). Each individual components is 
then convoluted with the electrostatic potential gener- 
ated by a sinusoidal density, as expressed in Eq. ([24|) 
to obtain the Fourier components v{Tj_;gz) of the open- 
boundary potential (step 2): 



v{rr.O) = -^J ln|r^-rY|p(rl;0)drl, 
v{r±;gz) = 2^ Ko{gz\r±-r^\)p{r^;gz)dr^ for g^ 0. 



(25) 
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Finally, the open-boundary potential is transformed back 
to real space (step 3). We underscore that this procedure 
directly extends the linear- average approximation since 
the linear average of the charge density corresponds to 
the first term of the one-dimensional Fourier decomposi- 
tion. Thus, averaging the charge density along the axis 
of periodicity amounts to restricting the Fourier series to 
its gz = term. 

To estimate errors resulting from this truncation, we 
analyze the asymptotic behavior of v{r±; 7^ 0) at large 
gz\r±\: 

v{Yr,gz) ^ J-—=== when gz\Y±_\ > 1. (26) 
V ^ V^^|r±| 

From Eq. [26l the validity of the linear average approach 
can be assessed by calculating the ratio of the cell size 
in the transverse direction L±_ (that is, the distance be- 
tween periodic replicas) to the typical wavelength Ay 
characterizing longitudinal inhomogeneities in the sys- 



tem. For large values of the dimensionless parameter 
periodic-image interactions are predominantly 
due to the logarithmic first-order contribution v{y^]{)) 
corresponding to the linear average of the charge density. 
Thus, as expected intuitively, the linear- average approx- 
imation is valid in this situation. In contrast, when Ay is 
comparable to the distance L±_ between periodic images, 
higher-order Fourier components v(r^; gz) corresponding 
to gz ~ 27r/A|| must also be taken into consideration. 

Despite its merit in discussing the validity of the linear- 
average approximation, determining the open-boundary 
potential using the preceding approach requires expen- 
sive summations for each point of the two-dimensional 
grid and for each longitudinal wavevector gz. Along the 
same methodological lines as those of the DCC algo- 
rithm, a substantial reduction of computational cost can 
be achieved by exploiting the periodic potential v' ^ whose 
longitudinal Fourier components can be computed inex- 
pensively using FFT techniques: 



Air 

v'{rr,9,) = ^^^p(g^+5.z)e^«-'-- for 5, ^ 0. 
tr^ SI + flS 



(27) 



After coarse-grid interpolation, the component of the 
open-boundary potential v{T_\_;gz) can be calculated at 
the boundaries of the domain, yielding Dirichlet bound- 
ary conditions for the smooth corrective components 
v^^^^{t_\_; gz) = v'{T_\_;gz) — v\T±;gz). The corresponding 
^2-dependent electrostatic problems read: 

r V2v^^^^(rx;0) = -47r(p) 

{ (28) 

[ (V^ - ^,2)^--(rx;^.) = for ^ 

These differential equations can be solved using efficient 
multigrid techniques. Once calculated, the longitudinal 
Fourier components of the electrostatic correction are 
added to those of the periodic potential, thereby recov- 
ering v{T_\_;gz). Finally, the potential v{t) is computed 
via an inverse Fourier transform. 



C. DCC Scheme for Two-dimensional Periodicity 

The electrostatic potential of a slab can be calculated 
in real space using a scheme similar to that presented 
above. The formalism is to a great extent analogous 
to that developed by Lang and Kohn for studying in- 
teractions between localized external charges and metal- 
lic surfaces [IT], and to the Green's function approach 



recently proposed by Otani and Sugino [42[. The pre- 
scription consists of performing two-dimensional Fourier 
transforms to obtain the charge-density profile p{z]g\\) 
associated with each wavevector gy = {gx^gy) parallel to 
the surface. Solving the electrostatic problem for sinu- 
soidal density layers, the two-dimensional Green's func- 
tions 5(z;g||) can be written as: 

g{z;0) = -2tt\zI 

C?(z;g||) = 27r for 511 7^0. ^^^^ 

^11 

Hence, as in the one-dimensional case, the density- 
average approximation is valid provided that the geomet- 
rical parameter L^/\\\ is large — this criterion is identical 
to that derived by Natan, Kronik, and Shapira [19]. In 
addition, the above expressions allow one to determine 
the corrective potential of a two-dimensional system by 
integrating the differential equations: 

^^--(^;0) = -4^(p) 

(30) 

= for ^11^0 

Parenthetically, it is important to note that Eq. [30] can 
be solved analytically, taking into account the boundary 
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FIG. 13: Total energy as a function of transverse cell size 
for a — [CH2CF2]3 — [CF2CH2]3— polymer chain without cor- 
rection, corrected using the density- countercharge scheme 
with full Fourier decomposition (DCC), and by limiting the 
density-countercharge decomposition to the linear- average 
g = component (DCC/LA). 
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FIG. 15: Longitudinal density response coefficient x{9\\) — 
dn{g\\) /dv{g\\) as a function of transverse cell size for a 
graphene sheet without correction, and corrected using the 
density-countercharge scheme with full Fourier decomposition 
(DCC). 
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FIG. 14: Force on one of the fluorine atoms along a trans- 
verse lattice direction as a function of transverse cell size for 
a — [CH2CF2]3 — [CF2CH2]3— polymer chain without correc- 
tion, corrected using the density-countercharge scheme with 
full Fourier decomposing (DCC), and by limiting the density- 
countercharge decomposition to the linear- average g = com- 
ponent (DCC/LA). 



conditions calculated by superposition — that is, by con- 
voluting the longitudinal components of Q and p (simi- 
larly to Eq. [25]), then subtracting out the components of 
v' . Therefore, the additional cost of the two-dimensional 
DCC correction is negligible. 



D. Applications 



The convergence of the total energy with respect 
to transverse cell size for a fluoropolymer chain 
-[CHsCFsls - [CF2CH2]3- of long periodicity Ay ^ 24 
bohr is depicted in Figure [131 We employ ultrasoft pseu- 
dopotentials [43] with energy cutoffs of 50 and 500 Ry 
for the plane- wave expansions of the electronic wavefunc- 



tions and charge density, respectively. The energy cutoff 
for calculating the corrective potential is 80 Ry. We use 
a shifted 1x1x2 mesh with cold-smearing occupations 
^6] (smearing temperature of 0.02 Ry). Within the lin- 
ear average approximation (DCC/LA), the corrected en- 
ergy closely coincides with the uncorrected energy due to 
the absence of polarization in the longitudinal average of 
the charge density. For the cell parameters considered, 
the geometrical ratio L^/\\\ varies from 0.5 to 0.9, that 
is, beyond the range of validity of the linear average ap- 
proximation. As a result, we observe that the DCC/LA 
energy converges slowly towards its asymptotic value. In 
contrast, the DCC scheme with full Fourier decomposi- 
tion significantly improves the convergence of the total 
energy (at 16 bohr, the accuracy of DCC energy is ap- 
proximately 5 X 10~^ Ry whereas that of the uncorrected 
and DCC/LA energies is approximately 10~^ Ry). Fig- 
ure [m depicts the convergence of the force on one of 
the fluorine atoms. Similarly to the convergence of the 
total energy, the atomic-force convergence is seen to im- 
prove substantially by applying the DCC correction: at 
16 bohr, the DCC force is converged within less than 
10~^ Ry/bohr, while that obtained without correction 
or using the DCC/LA scheme are converged within 10~^ 
Ry/bohr. We underscore that the additional computa- 
tional cost of the DCC correction is moderate. Indeed, 
at 16 bohr, the additional computational cost is ~8%. 

To conclude this study, we consider the electronic den- 
sity response of a graphene sheet subject to a perturba- 
tion field. Figure [T5l reports the dependence of the linear- 
response coefficient x{9\\) — dn{g\\) / dv{g\\) with respect 
to the interplane distance for a longitudinal sinusoidal 
perturbation of wavevector g\\ = ^ bohr~^. The wave- 
length of the perturbation field being large (Ay = 157 
bohr), the uncorrected response coefficient does not con- 
vergence until reaching cell sizes on the order of hun- 
dreds of bohrs. Contrary to uncorrected calculations, the 
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DCC-corrected linear response shows considerable con- 
vergence improvement with a negligible increase in com- 
putation cost. For comparison, at an interplane distance 
of = 50 bohr, the relative error in the uncorrected 
linear-response coefficient x{9\\) is on the order of 25%, 
while it is lower than 1% using the DCC correction. 
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V. CONCLUSION 

We have studied the analytical properties of the correc- 
tive potential, defined as the difference between the elec- 
trostatic potential and its periodic counterpart, unifying 
the Makov-Payne (FCC) and LMCC (GCC) schemes in 
the same class of periodic-image corrections and suggest- 
ing possible improvements for both methods. Based on 
these properties, we have shown that the periodic-image 
errors can be eliminated at a moderate computational 
cost of 0(M^/^), where M is the number of points of 
the mesh used in the calculation the corrective potential, 
which is generally about two orders of magnitude smaller 
than the number of points of the charge-density grid. 
The resulting density-countercharge (DCC) scheme owes 
its improved efficiency to the determination of the exact 
boundary conditions characterizing the electrostatic po- 
tential. In several cases of interest, we have shown that 
the DCC algorithm represents a beneficial compromise 
between cost and accuracy. The validity of the linear- 
and planar- average approximations routinely employed 
in the study of partially periodic systems has also been 
discussed. An efficient scheme going beyond these con- 
ventional approximations for inhomogeneous systems has 
been proposed and validated. 

Relevant applications for the DCC algorithm include 
the study of molecular adsorption at solid- vacuum inter- 
faces in the constant-charge regime, the determination of 
structural parameters, the correction of vibrational spec- 
tra, the inexpensive calculation of work functions, and 
the determination of linear-response properties with a 
reduced computational effort. 
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APPENDIX A: MADELUNG CONSTANTS AND 
GAUSSIAN POTENTIALS 

In this appendix, we determine the Madelung con- 
stants of periodic point charges immersed in a compen- 
sating jellium background in one, two, and three dimen- 
sions for lattices characterized by a single geometric pa- 
rameter L. A compilation of high-precision values for 
these fundamental constants is generally not found in the 
literature. 

These values are computed using the asymptotic ex- 
pansion of the Madelung constant a^r/L of an array of 
Gaussian charges of spread a in a compensating jellium, 
which is defined as: 



a./L = K(0)-<L(0))L'*-2, 



(Al) 



where d is the spatial dimension. To obtain the expansion 
of Q^cr/L ill the limit a / L 0, we may write v'^ ^(0) as: 



Wd 



(A2) 



(A3) 



where Vtd is the volume of d-dimensional unit cell, and 
= Lg denotes the dimensionless wave vector. Differen- 
tiating Wd with respect to jL^ ^ we obtain: 



= -^2^exp(-— ) 

/2 2 

= ^-^^exp(-^-^). (A4) 



In the limit a/L ^ 0, this derivative becomes: 



dwd 



d{ayL^) 



/ 2 \ n 



Integrating this expression, we obtain the asymptotic ex- 
pansions of '^^^(0) and a^r/L listed in Table [D 

Hence, the Madelung constant ao can be calculated 
with high accuracy from the expansion of a^r/L- In the 
case of a cubic lattice of point charges, we obtain: 
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FIG. 16: Convergence of the Madelung constant as a function 
of the geometric parameter L/cr for a cubic unit cell using the 
approximation given by Eq. IA6I (The black curve is Eq. IA6I 
without the ttct^/L^ and the complementary-error-function 
terms, the red curve is Eq. IA6I without the complementary- 
error- function term, and the blue curve is Eq. IA6|) . Note the 
negligible contribution of the complementary-error-function 
term beyond L/cr = 3 and the improvement in convergence 
brought about by the term ttct^/L^. 
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where n = (i, j, /c) denotes an integer vector. Figure [T6lil- 
lustrates the rapid convergence of the Madelung constant 
calculated from Eq. IA6I for a cubic cell. This expres- 
sion converges considerably faster than the expression 
frequently found in the literature: 



1 ^ 47r 



(A7) 



117^0 



Although a similar procedure can be applied without ad- 
ditional difficulty for any dimensionality, we draw atten- 
tion to the fact that in two dimensions, OLajL is not equal 
to the Madelung constant a in the limit a/L ^ 0, due to 



the logarithmic divergence of the potential. For a more 
complete discussion of the two-dimensional case, we re- 
fer the reader to the study of Cichocki and Felderhof 
[27]. As a final remark, we note that the one-dimensional 
Madelung constant can be determined analytically from 
the relation: 



+00 ^ 

^ TV' 

n=l 
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(A8) 
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stands for the Riemann zeta function. 
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FIG. 17: Accuracy of the total energy of a pyridazine cation 
as a function of computational time using the PCC, GCC, 
and MCC schemes for cell sizes in the range 12-19 bohr. The 
labels D (dipole) and Q (quadrupole) indicate the order of the 
multipole expansion. For each scheme the corrective potential 
is updated every five SCF iterations. 



APPENDIX B: PERFORMANCE OF THE 
MULTIPOLE-EXPANSION METHOD 



The performance of the multipole-expansion adapta- 
tion of the DCC scheme — the multipole-countercharge 
(MCC) correction — for a pyridazine cation is compared 
to that of the PCC and GCC schemes in Figure [TTl The 
size of the calculation cell ranges from 12 to 19 bohr. The 
parameters used in these calculations are those detailed 
in Sec. IIIIBI Note the good performance of the MCC 
approach, which improves the energy accuracy by almost 
one order of magnitude in comparison with the PCC and 
GCC schemes for cell sizes above 17 bohr. 
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